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Abstract 

Aggregation of particles whose interaction potential depends on their mutual orientation is con- 
sidered. The aggregation dynamics is derived using a version of Darcy's law and a variational 
principle depending on the geometric nature of the physical quantities. The evolution equations 
that result separate into two classes: either characteristic equations, or gradient flow equations. 
We derive analytical solutions of both types of equations which are collapsed (clumped) states and 
show their dynamical emergence from smooth initial conditions in numerical simuations. 
Keywords: gradient flows, blow-up, chemotaxis, parabolic-elliptic system, singular solutions 
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Introduction Many physical processes can be understood as aggregation of individual 
'components' at a variety of scales into a final 'product'. Diverse examples of such processes 
include the formation of stars, galaxies and solar systems at large scales, organization of 
insects and organisms into colonies at mesoscales and self-assembly of proteins, nanotubes 
or micro/nanodevices at micro- and nanoscales Some of these processes, such as nano- 
scale self-assembly of molecules are of great technological interest P, 0]. Of special interest 
is the case when the energy of interaction among the assembling particles depends not just 
on the distance, but also on mutual orientation. In particular, this Letter is motivated 
by recent experiments and simulations on self-assembly of non-circular particles (squares, 
hexagons etc.) J, 0, Q|. Due to the large number of particles involved in self-assembly 
(10 9 — 10 12 ), the development of continuum descriptions for aggregation or self-assembly is 
a natural approach toward theoretical understanding and modeling. 

Self-assembly at microscopic scales may be simpler than fluid dynamics, because the 
Reynolds number is so low that inertia is negligible. The classic examples of continuous 
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equations for aggregation are those of Debye-Hiickel p| and Keller-Segel (KS) jSj. For a 
recent review of developments in this area with an emphasis on biophysical modeling, see, 
for example Q. The physics of these models consists of a conservation law: d t p + divpu = 0, 
coupled with an evolution equation for velocity u which depends on the density p through a 
free energy E as u ~ pS75E/5p (velocity proportional to force), in which 'mobility' \i may 
so depend on the density. This relation is known as "Darcy's Law." In most recent works 
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and references within) the energy E is computed under the assumption that the 
interaction between the particles is central, i.e., that it does not depend on their orientation. 
This framework is simple and attractive. However, the self-assembly of some physical sys- 
tems may depend on geometric properties such as the mutual orientation of pieces. Examples 
of such systems range from micro-biological (mutual attraction of cells, viruses or proteins), 
to electromagnetic (dipoles in continuous media, orientation of domains), to interactions of 
living organisms (swarms, herds, flocks, etc.) 

The goal of this Letter is to formulate a procedure to derive and analyze evolution equations 
for systems that self assemble under a flow in which generalized 'velocity' is proportional to 
generalized 'force.' Although we aim to derive equations for orientation-dependent motion, 
we will also show that the procedure for obtaining such evolution equations applies generally 
for any type of continuum physical quantity. Thus, we obtain a family of evolution equations 
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for physical quantities such as active scalars, momenta and fluxes. (A complete list of equa- 
tions derived by our method is given in After these equations have been formulated, 
we discuss their most remarkable feature. Namely, they possess emergent solutions - or 
coherent structures - in which the self-assembly is singular (localized into delta functions). 
In addition, these localized singular solutions dominate the long-term dynamics. From the 
physical point of view, such localized, or quenched solutions would form the core of the pro- 
cesses of self-assembly and are therefore of great practical interest. Moreover, the formation 
of these localized solutions is driven by a combination of nonlinearity and nonlocality in an 
evolutionary process that admits an analytical description of their long-term dynamics. 
Geometric evolution equations The general principle for deriving an evolution equation 
for the self-assembly of any continuum physical quantity k may be stated as, "The local value 
of k remains invariant along the characteristic curves of a flow, whose velocity depends on 
k through an appropriate Darcy Law." This principle may be formulated in symbols as, 
dn(x(t),t)/dt = along dx/dt = u[k]. In the case of particle density p in n-dimensional 
space, for example, the number of particles n = p d n x has physical meaning. Hence, the time 
derivative of k in this case invokes the fundamental chain rule for the product of the density 
function times the volume element, k = p(x(t), t) d n x(t). Preservation of this product along 
dx/dt = u[p] yields 



in which the dependence of the velocity vector u[p] as a functional of p is yet to be de- 
termined. As mentioned above, the Darcy Law approach assumes that velocity u depends 
on density p through the gradient of the variation of free energy E (velocity proportional 
to 'force' with mobility /i). This assumption leads to the expected continuity equation for 
density, 



where sharp ( ■ )" denotes raising the vector index from covariant to contravariant, so its 
divergence may be taken. This gradient-flow method generalizes for other physical quantities 
with different geometrical meaning, by noticing that the invariance property dn(x(t), t)/dt = 
along dx/dt = u[k] may be expressed mathematically as dn/dt + £ U \ K ]K = 0, where £ u k 
is the Lie derivative with respect to the vector field u = u ■ V of any geometrical quantity 
The key question for the physical modeling is to identify the analog of Darcy's Law. 



(dtp + u[p] ■ Vp + p div u[p]) d n x(t) = , 



(1) 



d t p = - div p(pV5E/5p) s , 



(2) 
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Namely, what is the vector field u[k] when k is an arbitrary geometrical quantity? 
A surprising clue leading to the geometrical analog of Darcy's Law emerges when considering 
the spontaneous appearance of singularities in solutions of (E) for which p and SE/Sp depend 
on the average density, rather than its pointwise value Q]. Those singular solutions 
obey a weak form of the continuity equation (j2J), expressed by pairing it with an arbitrary 
smooth test function and integrating twice by parts, as 

(dtp , 0) = (5p , 5E/5p) to find (3) 
5p = — dw ppV(f) = — £ u (<i>) P so u(4>) = (pV4>f 

where spatial integration defines the real-valued pairing ( ■ , • ) between densities and their 
dual space of scalar functions. This expression for the weak solutions of the continuity 
equation (J2J) provides the clue we seek for expressing Darcy's Law for the self-assembly of 
an arbitrary geometric quantity k, not just a density. 

The clue we seek emerges upon re-doing the previous integrations by parts slowly and 
carefully as 

5E \ I 5E 

divp(/iV0)», — ) =: -/(/io0)», po — 

Here, we have introduced the diamond (o) as the dual of the Lie derivative under integration 
by parts for any dual pair (k, b) and any vector field u, (nob, u) = (k, —£ u b). This suggests 
u(0) = (pV(f>y — — (p O 0)". Using the definition of diamond again gives 

< | ,*) = ~ (0* o ^ > P o ^> = " <^f ^ , 0> (4) 

Replacing p k here generalizes (0) to any quantity k in an arbitrary vector space. We 
shall derive the particular versions of this equation for several physical quantities k. The 
result is a set of novel nonlocal characteristic equations which possess localized solutions. 
Next we shall outline a general procedure for finding those solutions. 

Singular solutions for the Geometric order parameter (GOP) equation (jl]). Nu- 
merical simulations show the development of singular solutions in k for some types of k 
(e.g., scalars) and the absence of those singularities for other k's (e.g., vector fields). These 



solutions are reminiscent of the clumpon singularities 
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3, 



which dominate the long term 
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dynamics once they have formed. Thus, we seek a particular solution of (JH) expressed as a 
5-function parametrized by coordinate(s) s on a submanifold of ambient space, 

n(x,t) = J p(s,t)5 {x — q(s,t)) ds . (5) 

To derive the equations for p(s,t) and q(s, t), we substitute (JHJ) into (JU) and integrate the 
right-hand side by parts to extract the term proportional to k as follows: 

|V(qM)) + ^-V<MqM))= UkA (6) 

where M K is a linear operator acting on <j) depending on the nature of k. For example, for 
densities k = pd 3 x, N K <t> = —jjiW5E/5p ■ V0(q(s,t)). If the right-hand side of © only 
contains the function <fi and its gradient, then singular solutions (0) are possible. However, 
one must be careful because this condition could be over-determined and thus the existence 
of singular solutions of (jlj) for the evolution of an arbitrary geometric quantity might not 
be guaranteed. Two classes of geometric quantities admitting singular solutions of (jlj) are 
known . One class includes, for example, scalars, 1-forms and 2- forms, and gives charac- 
teristic equations for which the characteristic velocity is a nonlocal vector function. A second 
class, which encompasses in particular an orientation-dependent density, is a nonlinear non- 
local diffusion equation. In each case, one must compute the particular expression for the 
right hand side of (jlj), and integrate by parts to extract the function <fi and its derivatives. 
Since our main interest lies with the second case, we shall describe the nonlocal character- 
istic equations only very briefly. It is essential to discuss them, however, because of their 
interesting mathematical and physical properties and their many possible applications. 
Nonlocal characteristic equations The fundamental example is for a scalar k = f. The 
evolution of a scalar by (jlj) obeys 

f = -^Jf).M/] = -(fvM/])".V/. (7) 

Equation (J7J) can be rewritten in characteristic form as df /dt = on d^jdt = ^|jV/i[/] 
The characteristic speeds of this equation are nonlocal when SE/5f and fi are chosen to 
depend on the average value, f. 

One may verify that the scalar equation (j?J) admits weak solutions Figure 1 shows the 
spatio-temporal numerical evolution of H * f given by (J7J) with initial conditions of the type 
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dPa(t,s) , . /SE 

— — — = p a (t, s) div (— V/i[/j 
dq (t,g) /<5£ \t 

Pa(M) ^ = p a (t,s) V/i[/]J 



(J3J) with ^-functions whose strengths are random numbers between ±1/8. We have taken 
SE/Sf = H * / where H is the inverse Helmholtz operator H(x) = e~^ x ^ a with a — 1. 
We see the evolution of sharp ridges in / = H * /, which corresponds to 5-functions in the 
solutions f(x, t). 

Explicit equations for the evolution of strengths p a and coordinates q a for a sum of S- 
functions in (jSJ) may be derived using (JBJ) as 

(8) 

7 a==g (t,s) 

for a = 1, 2, . . . , N. A solution containing a single 5-function satisfies p = — Ap 3 , so an 
initial condition p(0) = po, evolves according to l/p(t) 2 = 1/pq + 2 At. For the choice of 
our parameters in simulations, A = 2. The comparison of 1/p 2 from numerics with this 
theoretical prediction is shown in Figure 2. 

We only mention that the evolution of a general 1-form (velocity) and 2-form (flux) may also 
be cast into the form of a nonlocal characteristic equation. For a 1-form in two- or three- 
dimensional space A • dx, the equation is quite complicated However, that equation 
allows a simplification in the case A = V^. The equation for potential ip reads, with 
nonlocal SE/Stp and 

£-(»v*i)'.v». do, 

This equation for the potential i/j has the same nonlocal characteristic structure as the scalar 
equation J7J). 

Similarly, the evolution of 2-form fluxes B • dS = Bidx 2 A dx 3 — i^dx 1 A dx 3 + f^dx 1 A dx 2 
also simplifies. Again, the evolution for a general flux B • dS is quite complicated, but it can 
be simplified for the case divB = when B only depends on two coordinates (x,y). In this 
case, equation (jlj) may be written for the stream function \I/ in B = curl = x z as 

5* _ /SE 



at v» v *J - v *> < n > 

when we choose mobility to be \x = curl (z$) = V$ x z. Choosing SE/5^ and $ to depend 
on the average value ^ again yields a nonlocal characteristic equation. 
Nonlocal aggregation equations Next, we consider the case which is the main focus of 
the paper, namely the orientation-dependent density. The evolution equations arising in 
this case are similar to (0) for mass density. 
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Density at each point is given by one number, whereas the number of coordinates necessary to 
describe the orientation by an element of the rotation group SO(n) depends on the dimension 
of the space n. In two dimension, one number - turning angle - is enough, whereas in three 
dimensions three numbers are required. Keeping in spirit with the framework (@J), we see 
that we cannot make a theory describing evolution of orientation per se, because SO(n) 
is not a linear space. Instead, we need to consider the tangent space to SO(n) taken at 
the identity, commonly denoted so(n), which is a linear space. The description in terms 
of coordinates on so{n) is simple and corresponds to physical intuition. So, let us define 
n to be a pair of densities (p, a) where p is mass or number density and a is orientation 
density of particles, taking values in so(n). Hence, we set k = ^2K a e a where e a are basis 
vectors spanning the space of k, a — 0, 1, 2, 3 with a = corresponding to the mass density. 
In 3-dimensional space, an elegant description of k exists using quaternions, obtained by 
associating the density part a = with the scalar part and a = 1, 2, 3 with the 'vector' part 
of the quaternion. In two dimensions, an analogous description uses complex numbers. For 
orient at ion- dependent aggregation equation in 3D equation (JH) reads 

^ = -div(/[K]( K a V^)), a, b = 0,1, 2, 3, (12) 

or, explicitly for p and a b , b = 1,2, 3: 

(Summation on a is assumed, but no summation on b.) 

Interestingly enough, equations (|13ll4j) allow weak solutions of the form (JHJ). A special 
condition for existence of such solutions is that all mobilities p b must be constant along 
the flow. Then, even though ansatz (jSJ) leads to an over-determined system, this system is 
consistent and for any real K, the following is an exact solution of (141 

p „ ,M l_« 8 SE 



A'¥=(f>"## V- (15) 

r L 1 dt \ dxi 5n a r=o(s,t)J v ' 



dt V dxi 8n a r=q(s,t); 
In Figure 3, we present a simulation of the density and orientation in two dimensions. 
Starting with density and orientation spread unequally over two Gaussian 'clumps', we end 
up with a stationary solution whose density and orientation are concentrated into one clump. 
This is the nature of self assembly in this case. 



Further developments Equation (j3J) also encompasses ideal fluid vorticity dynamics 
(which is not associated with Darcy's law). For the pairing (oj , m) — j uj • mdV be- 
tween a divergenceless vector field uj = uj ■ V and a one-form density m = m • dx ® dV, 
set i raj = (Su , ^\ and choose variation <5cj = «£ cur i m k; = [curlm , a;]. After two 
integrations by parts, one finds d t to = — [ curl^j , uj ] which recovers fluid vortex dynamics 
when curl|^ = u and uj = curlu. It will be interesting to see what other physical systems 
remain to be discovered in equation (j3J). 
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Figure Captions. 



Figure 1. Numerical simulation of scalar /, starting with initial condition which is a set of 5- 
functions in /. The vertical coordinate represents f = H*f, which remains finite even 
when / forms ^-functions, where H(x) = e~' x L The horizontal coordinate is space; 
the vertical coordinate is time. 

Figure 2. Evolution of the 5-function strength l/p(t) 2 versus time (circles). The theoretical 
prediction + At is shown as a solid line obtained without any fitting parameters. 

Figure 3. Numerical simulation of density p starting with initial conditions for p and a that are 

2 ill 

two Gaussian clumps of unequal strength (p, <r)(r, 0) = Y2dumps(Poi a o) e F (p<T) - The 
initial conditions for p and a are similar in shape but different in width. For this 
simulation, we have taken l a = 21 p , H(pc) = e~' x L Top: initial conditions, bottom: 
finite solution for t = 30. Left column: density, right column: orientation. White 
represents the domains of high density/orientation. 
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